//
//  Move.h
//  phi_four
//
//  Created by Snir Gazit on 3/9/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#ifndef Move_h
#define Move_h
#include "../HelpLibs/RandomGen.h"
#include "Lattice.h"
#include "../HelpLibs/QmcParameters.h"
namespace ba = boost::accumulators;

typedef RandomVariable<MTGen, UniReal> UniRealVariable;
typedef RandomVariable<MTGen, UniInt> UniIntVariable;
typedef RandomVariable<MTGen, Binomial> BinomialVariable;
typedef RandomVariable<MTGen, Discrete> DiscreteVariable;
typedef UniRealVariable::param_type UniRealVariableParams;
typedef UniIntVariable::param_type UniIntVariableParams;
typedef BinomialVariable::param_type BinomialVariableParams;
typedef DiscreteVariable::param_type DiscreteVariableParams;

using namespace std;

class Move {
protected:
	MTGen *eng_;
	Lattice *lat_;
	UniRealVariable zeroOne_;
public:
	void init(MTGen * eng, Lattice *lat);
	bool makeMove();
	bool accept(double p);
	bool checkMove();
};

// Worm Alg Moves
template<class BondType>
class ShiftMove: public Move {
	typedef boost::array<BondType*, 2 * DIM> BondArray;
	UniIntVariable chooseMove_;
	BoltzWeightType CreateRatio_, RemoveRatio_;
	double g_;
public:
	void init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
			BoltzWeightType boltzRatio);
	bool makeMove();
	int& getBondVal(BondArray &bondArray, int dir, CreateRemove createRemove);
	void setCoupeling(BoltzWeightType boltzRatio, double g);
};

template<class BondType>
class CreatePairMove: public Move {
	typedef boost::array<BondType*, 2 * DIM> BondArray;
	UniIntVariable chooseSite_;
	UniIntVariable chooseMove_;
	BoltzWeightType boltzRatio_;
	double g_;
public:
	void init(MTGen *eng, Lattice *lat, QmcParameters &qmcParameters,
			BoltzWeightType boltzRatio);
	void setCoupeling(BoltzWeightType boltzRatio, double g);
	bool makeMove();
	int& getBondVal(BondArray &bondArray, int dir, CreateRemove createRemove);
};

// Shift Move

template<class BondType>
void ShiftMove<BondType>::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
	chooseMove_.init(eng, UniIntParams(0, 2 * 2 * DIM - 1));
//	createRemove_.init(eng, BinomialVariableParams(1, 0.5));
	//openCloseRatio_ = qmcParameters.openCloseRatio.getVal();
	setCoupeling(boltzRatio, qmcParameters.g.getVal());
}
template<class BondType>
void ShiftMove<BondType>::setCoupeling(BoltzWeightType boltzRatio, double g) {
	g_ = g * (gFactor<BondType>());
	RemoveRatio_[0] = 0.;
	for (int j = 0; j < MAX_BOND; j++) {
		CreateRatio_[j] = boltzRatio[j] * g_;
		if (j != (MAX_BOND - 1))
			RemoveRatio_[j + 1] = 1 / (boltzRatio[j] * g_);
	}
}
template<class BondType>
bool ShiftMove<BondType>::makeMove() {
// Get Ira
	Lattice::SiteIter ira = lat_->getIra();
	// Choose Move
	int move = chooseMove_.draw();
	// Create Or Remove
	CreateRemove createRemove = (move < 2 * DIM) ? CREATE : REMOVE;
	// Direction
	int dir = move % (2 * DIM);
	Lattice::SiteIter nSite = ira->nSiteArray_[dir];
	BondArray& bondArray = ira->getBondArray<BondType>();
	int &bondVal = getBondVal(bondArray, dir, createRemove);
	// Create Move
	if (createRemove == CREATE) {
		// Annihilate Worm
		if (nSite == lat_->getMasha()) {
			double p = g_ / double(bondVal + 1);
			if (this->accept(p)) {
				bondVal++;
				// lat_->setIra(nSite); (needed??)
				lat_->sectorType_ = ZSECTOR;
				return true;
			}
			// Create New Bond
		} else {
			// Get Total Bonds
			int &totalBondsType = nSite->getTotalBondsType<BondType>();
			int totalBondsOther = nSite->getTotalBondsOther<BondType>();
			BondArray bondArray = ira->getBondArray<BondType>();
			int &bondVal = getBondVal(bondArray, dir, createRemove);
			double p = (CreateRatio_[totalBondsType + totalBondsOther]
					* addRatio<BondType>(totalBondsType, totalBondsOther))
					/ double(bondVal + 1);
			if (this->accept(p)) {
				bondVal++;
				totalBondsType++;
				lat_->setIra(nSite);
				assert(this->checkMove());
				return true;
			}
		}
	} else if (bondVal != 0) {
		// Annihilate Worm
		if (nSite == lat_->getMasha()) {
			// Get Total Bonds Ira
			int &totalBondsTypeIra = ira->getTotalBondsType<BondType>();
			int totalBondsOtherIra = ira->getTotalBondsOther<BondType>();
			// Get Total Bonds Masha
			int &totalBondsTypeMasha = nSite->getTotalBondsType<BondType>();
			int totalBondsOtherMasha = nSite->getTotalBondsOther<BondType>();
			double removeIra = RemoveRatio_[totalBondsTypeIra
					+ totalBondsOtherIra]
					* removeRatio<BondType>(totalBondsTypeIra,
							totalBondsOtherIra);
			double removeMasha = RemoveRatio_[totalBondsTypeMasha
					+ totalBondsOtherMasha]
					* removeRatio<BondType>(totalBondsTypeMasha,
							totalBondsOtherMasha);
			double p = g_ * removeIra * removeMasha * double(bondVal);
			if (this->accept(p)) {
				bondVal--;
				assert(bondVal < MAX_BOND);
				totalBondsTypeIra--;
				totalBondsTypeMasha--;
				lat_->sectorType_ = ZSECTOR;
				// lat_->setIra(nSite); (Needed?)
				return true;
			}
		}
		// Remove Bond
		else {
			// Get Total Bonds
			int &totalBondsType = ira->getTotalBondsType<BondType>();
			int totalBondsOther = ira->getTotalBondsOther<BondType>();
			double p = (RemoveRatio_[totalBondsType + totalBondsOther]
					* removeRatio<BondType>(totalBondsType, totalBondsOther))
					* double(bondVal);
			if (this->accept(p)) {
				bondVal--;
				assert(bondVal < MAX_BOND);
				totalBondsType--;
				lat_->setIra(nSite);
				assert(this->checkMove());
				return true;
			}
		}
	}
	return false;
}
template<>
int& ShiftMove<XYBondType>::getBondVal(
		boost::array<XYBondType*, 2 * DIM> &bondArray, int dir,
		CreateRemove createRemove) {
	if (dir < DIM) {
		if (createRemove == CREATE)
			return (*(bondArray)[dir])[0];
		else
			return (*(bondArray)[dir])[1];
	} else {
		if (createRemove == CREATE)
			return (*(bondArray)[dir])[1];
		else
			return (*(bondArray)[dir])[0];

	}
}
template<>
int& ShiftMove<IsingBondType>::getBondVal(
		boost::array<IsingBondType*, 2 * DIM> &bondArray, int dir,
		CreateRemove createRemove) {
	return (*(bondArray)[dir]);
}
// Create Pair Move
template<class BondType>
void CreatePairMove<BondType>::init(MTGen *eng, Lattice *lat,
		QmcParameters &qmcParameters_, BoltzWeightType boltzRatio) {
	static_cast<Move *>(this)->init(eng, lat);
//boltzRatio_ = boltzRatio;
	chooseSite_.init(eng, UniIntParams(0, lat->getNumOfSites() - 1));
	chooseMove_.init(eng, UniIntParams(0, 2 * 2 * DIM - 1));
	boltzRatio_ = boltzRatio;
}

template<class BondType>
bool CreatePairMove<BondType>::makeMove() {
	// JumpSite
	int jumpSite = chooseSite_.draw();
	Lattice::SiteIter jump = lat_->getSite(jumpSite);
	// Move
	int move = chooseMove_.draw();
	// Create Or Remove
	CreateRemove createRemove = (move < 2 * DIM) ? CREATE : REMOVE;
	// Direction
	int dir = move % (2 * DIM);
	Lattice::SiteIter nSite = jump->nSiteArray_[dir];
	BondArray &bondArray = jump->getBondArray<BondType>();
	int &bondVal = getBondVal(bondArray, dir, createRemove);
	// Create Move
	if (createRemove == CREATE) {
		// Get Total Bonds
		// Get Total Bonds Ira
		int &totalBondsTypeMasha = jump->getTotalBondsType<BondType>();
		int totalBondsOtherMasha = jump->getTotalBondsOther<BondType>();
		// Get Total Bonds Masha
		int & totalBondsTypeIra = nSite->getTotalBondsType<BondType>();
		int totalBondsOtherIra = nSite->getTotalBondsOther<BondType>();
		double addIra = boltzRatio_[totalBondsTypeIra + totalBondsOtherIra]
				* addRatio<BondType>(totalBondsTypeIra, totalBondsOtherIra);
		double addMasha =
				boltzRatio_[totalBondsTypeMasha + totalBondsOtherMasha]
						* addRatio<BondType>(totalBondsTypeMasha,
								totalBondsOtherMasha);
		double p = g_ * addIra * addMasha / double(bondVal + 1);
		if (this->accept(p)) {
			bondVal++;
			totalBondsTypeIra++;
			totalBondsTypeMasha++;
			lat_->setIra(nSite);
			lat_->setMasha(jump);
			// lat_->setIra(nSite); (needed??)
			lat_->sectorType_ = getSectorType<BondType>();
			return true;
			assert(this->checkMove());
		}
	} else {
		double p = (1 / g_) * double(bondVal);
		if (this->accept(p)) {
			bondVal--;
			lat_->setIra(nSite);
			lat_->setMasha(jump);
			lat_->sectorType_ = getSectorType<BondType>();
			return true;
			assert(this->checkMove());
		}
	}
	return false;
}
template<class BondType>
void CreatePairMove<BondType>::setCoupeling(BoltzWeightType boltzRatio,
		double g) {
	g_ = g * (gFactor<BondType>());
	boltzRatio_ = boltzRatio;
}
template<>
int& CreatePairMove<XYBondType>::getBondVal(
		boost::array<XYBondType*, 2 * DIM> &bondArray, int dir,
		CreateRemove createRemove) {
	if (dir < DIM) {
		if (createRemove == CREATE)
			return (*(bondArray)[dir])[0];
		else
			return (*(bondArray)[dir])[1];
	} else {
		if (createRemove == CREATE)
			return (*(bondArray)[dir])[1];
		else
			return (*(bondArray)[dir])[0];

	}
}
template<>
int& CreatePairMove<IsingBondType>::getBondVal(
		boost::array<IsingBondType*, 2 * DIM> &bondArray, int dir,
		CreateRemove createRemove) {
	return (*(bondArray)[dir]);
}
#endif
